function [a0, a1, tau] = expfit(x, y)
% Perform a single exponential fit to data using Chebyshev polynomial method.
% Equation fit: y = a1 * exp(-x/tau) + a0
% Call: [a0 a1 tau] = expfit(x,y);
% Calling parameter x is the time base, y is the data to be fit.
% Returned values: a0 is the offset, a1 is the amplitude, tau is the time
% constant (scaled in units of x).
% Relies on routines chebftd to generate polynomial coeffs, and chebint to compute the
% coefficients for the integral of the data. These are now included in this
% .m file source.
% This version is based on the one in the pClamp manual: HOWEVER, since
% I use the bounded [-1 1] form for the Chebyshev polynomials, the coefficients are different,
% and the resulting equation for tau is different. I manually optimized the tau
% estimate based on fits to some simulated noisy data. (Its ok to use the whole range of d1 and d0 
% when the data is clean, but only the first few coeffs really hold the info when
% the data is noisy.)
% NOTE: The user is responsible for making sure that the passed data is appropriate,
% e.g., no large noise or electronic transients, and that the time constants in the
% data are adequately sampled.
% To do a double exp fit with this method is possible, but more complex.
% It would be computationally simpler to try breaking the data into two regions where
% the fast and slow components are dominant, and fit each separately; then use that to
% seed a non-linear fit (e.g., L-M) algorithm.
% Final working version 4/13/99 Paul B. Manis
%
n = 30; % default number of polynomials coeffs to use in fit
a = min(x);
b = max(x);
d0 = chebftd(a, b, n, x, y); % coeffs for data trace...
d1 = chebint(a, b, d0, n); % coeffs of integral...
%d2 = chebint(a, b, d1, n);
%d1./d0
tau =-mean(d1(3:4)./d0(3:4));
g = exp(-x/tau);
dg = chebftd(a, b, n, x, g); % generate chebyshev polynomial for unit exponential function
%% now estimate the amplitude from the ratios of the coeffs.
a1 = estimate(d0, dg, 2);
% get the offset here
a0 = (d0(1)-a1*dg(1))/2;
%fprintf(1, 'a0: %12.5f  a1: %12.5f   tau: %8.3f\n', a0, a1, tau);
return

function [a] = estimate(c, d, m)
% compute optimal estimate of parameter from arrays of data
n = length(c);
a = sum(c(m:n).*d(m:n))/sum(d(m:n).^2);
return;


% note that the following function has been translated into c code and exists as a DLL...
% 9/12/02 P. Manis
% to reutilize the present routine, change the routine name back to chebftd (which will override)..
%
function c = chebftdx(a, b, n, t, d)
% Chebyshev fit; from Press et al, p 192.
% matlab code P. Manis 21 Mar 1999
% "Given a function func, lower and upper limits of the interval [a,b], and
% a maximum degree, n, this routine computes the n coefficients c[1..n] such that
% func(x) sum(k=1, n) of ck*Tk(y) - c0/2, where y = (x -0.5*(b+a))/(0.5*(b-a))
% This routine is to be used with moderately large n (30-50) the array of c's is 
% subsequently truncated at the smaller value m such that cm and subsequent
% terms are negligible."
%
% This routine is modified so that we find close points in x (data array) - i.e., we find
% the best Chebyshev terms to describe the data as if it is an arbitrary function.
% t is the x data, d is the y data...
fprintf('\nEntering Chebyshev\n');
bma = 0.5*(b-a);
bpa = 0.5*(b+a);
inc = t(2)-t(1);
for k=0:n-1
   y = cos(pi*(k+0.5)/n);
   pos = round((y*bma+bpa)/inc);
   if(pos < 1) pos = 1;
   end
   if(pos >= length(d-1)) pos = length(d)-1;
   end
   %  disp(sprintf('Pos: %8d, equiv x: %8.3f', pos, y*bma+bpa))
   f(k+1) = d(pos+1);
%   fprintf('pos: %d d(pos+1): %f f(%d+1): %f\n', pos, d(pos+1), k, f(k+1));
end
fac = 2/n;
%fprintf('fac: %f bma=%f bpa=%f inc=%f n=%d\n', fac, bma, bpa, inc, n);
for j=0:n-1
   sum=0;
   for k=0:n-1
      sum = sum+ f(k+1)*cos(pi*j*(k+0.5)/n);
   end
   c(j+1)=fac*sum;
%	fprintf('j = %d, sum = %f, c[j] = %f\n', j+1, sum, c(j+1));
end

function [cint] = chebint(a, b, c, n)
% Given a, b, and c[1..n] as output from chebft or chebftd, and given n,
% the desired degree of approximation (length of c to be used),
% this routine computes cint, the Chebyshev coefficients of the
% integral of the function whose coeffs are in c. The constant of
% integration is set so that the integral vanishes at a.
% Coded from Press et al, 3/21/99 P. Manis
%
sum = 0;
fac = 1.;
con = 0.25*(b-a); % factor that normalizes the interval
for j=2:n-1
   cint(j)=con*(c(j-1)-c(j+1))/(j-1);
   sum = sum + fac * cint(j);
   fac = - fac;
end
cint(n) = con*c(n-1)/n;
sum = sum + fac*cint(n);
cint(1) = 2.0*sum; % set constant of integration.
return;
